%% pde_MMF_t
clc;
clear;

index = 1;

syms x y t nu nu_r S nu_m C0 Ca Cd
switch index
    case 1
        u1 = pi*sin(pi*x).^2.*sin(2*pi*y).*sin(t);
        u2 = -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        B1 = sin(pi*x).*cos(pi*y).*cos(t);
        B2 = -cos(pi*x).*sin(pi*y).*cos(t);
        w3 = pi*sin(pi*x).^2.*sin(pi*y).^2.*sin(t);
        p = 10*cos(pi*x).*cos(pi*y).*sin(t);
end

u1dx = diff(u1,x);
u1dy = diff(u1,y);
u1dt = diff(u1,t);
u2dx = diff(u2,x);
u2dy = diff(u2,y);
u2dt = diff(u2,t);
B1dx = diff(B1,x);
B1dy = diff(B1,y);
B1dt = diff(B1,t);
B2dx = diff(B2,x);
B2dy = diff(B2,y);
B2dt = diff(B2,t);
w3dx = diff(w3,x);
w3dy = diff(w3,y);
w3dt = diff(w3,t);
pdx = diff(p,x);
pdy = diff(p,y);
u1dxdx = diff(u1dx,x);
u1dydy = diff(u1dy,y);
u2dxdx = diff(u2dx,x);
u2dydy = diff(u2dy,y);
w3dxdx = diff(w3dx,x);
w3dydy = diff(w3dy,y);
B1dxdx = diff(B1dx,x);
B1dydy = diff(B1dy,y);
B1dxdy = diff(B1dx,y);
B1dydx = diff(B1dy,x);
B2dxdx = diff(B2dx,x);
B2dydy = diff(B2dy,y);
B2dxdy = diff(B2dx,y);
B2dydx = diff(B2dy,x);
u1B2dx = diff(u1*B2,x);
u1B2dy = diff(u1*B2,y);
u2B1dx = diff(u2*B1,x);
u2B1dy = diff(u2*B1,y);

f11 = u1dt - (nu+nu_r)*(u1dxdx+u1dydy) + (u1*u1dx+u2*u1dy) + pdx + S*B2*(B2dx-B1dy) - 2*nu_r*w3dy;
f12 = u2dt - (nu+nu_r)*(u2dxdx+u2dydy) + (u1*u2dx+u2*u2dy) + pdy + S*B1*(B1dy-B2dx) - 2*nu_r*(-w3dx);
f2 = w3dt -(Ca+Cd)*(w3dxdx+w3dydy) + (u1*w3dx+u2*w3dy) + 4*nu_r*w3 - 2*nu_r*(u2dx-u1dy);
f31 = B1dt + nu_m*(B2dxdy-B1dydy) - (u1B2dy-u2B1dy);
f32 = B2dt + nu_m*(B1dydx-B2dxdx) - (u2B1dx-u1B2dx);

disp("===u1导数===");
matlabFunction(simplify(u1dx), "Vars", {x,y,t})
matlabFunction(simplify(u1dy), "Vars", {x,y,t})
disp("===u2导数===");
matlabFunction(simplify(u2dx), "Vars", {x,y,t})
matlabFunction(simplify(u2dy), "Vars", {x,y,t})
disp("===B1导数===");
matlabFunction(simplify(B1dx), "Vars", {x,y,t})
matlabFunction(simplify(B1dy), "Vars", {x,y,t})
disp("===B2导数===");
matlabFunction(simplify(B2dx), "Vars", {x,y,t})
matlabFunction(simplify(B2dy), "Vars", {x,y,t})
disp("===w3导数===");
matlabFunction(simplify(w3dx), "Vars", {x,y,t})
matlabFunction(simplify(w3dy), "Vars", {x,y,t})
disp("===p导数===");
matlabFunction(simplify(pdx), "Vars", {x,y,t})
matlabFunction(simplify(pdy), "Vars", {x,y,t})
disp("===非定常MMF方程===")
matlabFunction(simplify(f11), "Vars", {x y t nu nu_r S nu_m C0 Ca Cd})
matlabFunction(simplify(f12), "Vars", {x y t nu nu_r S nu_m C0 Ca Cd})
matlabFunction(simplify(f2), "Vars", {x y t nu nu_r S nu_m C0 Ca Cd})
matlabFunction(simplify(f31), "Vars", {x y t nu nu_r S nu_m C0 Ca Cd})
matlabFunction(simplify(f32), "Vars", {x y t nu nu_r S nu_m C0 Ca Cd})
